Setup

Load R libraries

library(data.table)
library(ggplot2)
library(ggpubr)
library(tidyr)
library(limma)
library(stringr)

theme_set(theme_classic())

cell_type_name = params$cell_type_name
graph_weight = params$graph_weight

cell_type_name
## [1] "cd8"
graph_weight
## [1] "1.0"

Check enrichment of gene sets

Read in gene info and gene set assignments

file_tag = sprintf("%s_BL_%s", cell_type_name, graph_weight)

assayed_genes = scan(sprintf("output/gene_list_%s.txt", file_tag), 
                     what = character(), sep="\n")

gene_sets = scan(sprintf("output/name_s_%s.txt", file_tag), 
                 what = character(), sep="\n")

gene_sets = sapply(gene_sets, strsplit, USE.NAMES=FALSE, split=",")
n_genes   = sapply(gene_sets, length)
names(n_genes) = NULL
summary(n_genes)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   14.00   15.00   14.35   15.00   17.00
length(n_genes)
## [1] 40
sort(n_genes)
##  [1]  5 11 12 13 13 13 14 14 14 14 14 14 14 14 14 14 14 14 15 15 15 15 15 15 15
## [26] 15 15 15 15 15 15 15 15 15 16 16 16 17 17 17

Find gene symbols using the alias2Symbol function from limma.

a2s = rep(NA, length(assayed_genes))
for(i in 1:length(assayed_genes)){
  gi = assayed_genes[i]
  ai = alias2Symbol(gi)
  if(length(ai) > 1){
    print(gi)
    print(ai)
  }
  a2s[i] = ai[1]
}
## [1] "QARS"
## [1] "EPRS1" "QARS1"
## [1] "SEPT2"
## [1] "SEPTIN6" "SEPTIN2"
table(is.na(a2s))
## 
## FALSE  TRUE 
##  1607    42
table(a2s == assayed_genes, useNA = 'ifany')
## 
## FALSE  TRUE  <NA> 
##    42  1565    42
gene_info = data.table(sym_in_data = assayed_genes, sym_limma = a2s)

gene_info[sym_in_data != sym_limma,]
##     sym_in_data sym_limma
##  1:    C10orf91 LINC02870
##  2:    C12orf10      MYG1
##  3:    C12orf45  NOPCHAP1
##  4:     C6orf48    SNHG32
##  5:     C6orf99 LINC02901
##  6:    CXorf40A     EOLA1
##  7:     CXorf57      RADX
##  8:     FAM102A     EEIG1
##  9:     FAM173A    ANTKMT
## 10:     FAM213B    PRXL2B
## 11:       H2AFX      H2AX
## 12:   HIST1H2AG    H2AC11
## 13:   HIST1H2BK    H2BC12
## 14:   HIST1H2BN    H2BC15
## 15:    HIST1H3A      H3C1
## 16:    HIST1H3H     H3C10
## 17:    HIST1H4C      H4C3
## 18:   HIST2H2BF    H2BC18
## 19:    KIAA0391     PRORP
## 20:        QARS     EPRS1
## 21:       SEPT6   SEPTIN6
## 22:       ARNTL     BMAL1
## 23:    C12orf65     MTRFR
## 24:    C16orf72   HAPSTR1
## 25:      CCDC84   CENATAC
## 26:      DOPEY2     DOP1B
## 27:     FAM126B     HYCC2
## 28:    FAM160B1    FHIP2A
## 29:        H1FX     H1-10
## 30:       H2AFJ      H2AJ
## 31:       HEXDC      HEXD
## 32:    HIST1H1C      H1-2
## 33:    HIST1H1D      H1-3
## 34:    HIST1H1E      H1-4
## 35:    KIAA1109     BLTP1
## 36:    KIAA1551     RESF1
## 37:        MKL1     MRTFA
## 38:       NARFL     CIAO3
## 39:       SEPT2   SEPTIN6
## 40:      TARSL2     TARS3
## 41:      TMEM8A     PGAP6
## 42:       WDR60   DYNC2I1
##     sym_in_data sym_limma
gene_info[, gene_symbol := sym_in_data]
gene_info[which(sym_in_data != sym_limma & (gene_symbol != "MT-CO2")), 
                gene_symbol := sym_limma]

dim(gene_info)
## [1] 1649    3
gene_info[1:5,]
##    sym_in_data sym_limma gene_symbol
## 1:      ABLIM1    ABLIM1      ABLIM1
## 2:  AC004687.1      <NA>  AC004687.1
## 3:  AC004854.2      <NA>  AC004854.2
## 4:  AC007384.1      <NA>  AC007384.1
## 5:  AC007952.4      <NA>  AC007952.4
t1 = table(gene_info$gene_symbol)
table(t1)
## t1
##    1    2 
## 1647    1
gene_info[gene_symbol %in% names(t1)[t1 == 2],]
##    sym_in_data sym_limma gene_symbol
## 1:       SEPT6   SEPTIN6     SEPTIN6
## 2:       SEPT2   SEPTIN6     SEPTIN6
gene_info[sym_in_data == "HIST1H2BC", gene_symbol:="H2BC4"]
gene_info[sym_in_data == "HIST1H2BG", gene_symbol:="H2BC8"]
gene_info[sym_in_data == "SEPT6", gene_symbol:="SEPTIN6"]
gene_info[sym_in_data == "SEPT2", gene_symbol:="SEPTIN2"]

Read in cell type-specific expression data

Cell type specific gene expression were downloaded from protein atlas.

ct = fread("../Annotation/rna_single_cell_type.tsv.gz")
dim(ct)
## [1] 1587110       4
ct[1:5,]
##               Gene Gene name             Cell type  nTPM
## 1: ENSG00000000003    TSPAN6            Adipocytes 149.5
## 2: ENSG00000000003    TSPAN6 Alveolar cells type 1   6.1
## 3: ENSG00000000003    TSPAN6 Alveolar cells type 2  10.5
## 4: ENSG00000000003    TSPAN6            Astrocytes  13.9
## 5: ENSG00000000003    TSPAN6               B-cells   1.5
length(unique(ct$`Cell type`))
## [1] 79
table(ct$`Cell type`)
## 
##                      Adipocytes           Alveolar cells type 1 
##                           20090                           20090 
##           Alveolar cells type 2                      Astrocytes 
##                           20090                           20090 
##                         B-cells             Basal keratinocytes 
##                           20090                           20090 
##           Basal prostatic cells         Basal respiratory cells 
##                           20090                           20090 
## Basal squamous epithelial cells                   Bipolar cells 
##                           20090                           20090 
##          Breast glandular cells      Breast myoepithelial cells 
##                           20090                           20090 
##                  Cardiomyocytes                  Cholangiocytes 
##                           20090                           20090 
##                      Club cells           Collecting duct cells 
##                           20090                           20090 
##        Cone photoreceptor cells                Cytotrophoblasts 
##                           20090                           20090 
##                 dendritic cells              Distal enterocytes 
##                           20090                           20090 
##            Distal tubular cells                    Ductal cells 
##                           20090                           20090 
##                Early spermatids      Endometrial ciliated cells 
##                           20090                           20090 
##       Endometrial stromal cells               Endothelial cells 
##                           20090                           20090 
##           Enteroendocrine cells                 Erythroid cells 
##                           20090                           20090 
##              Excitatory neurons        Exocrine glandular cells 
##                           20090                           20090 
##       Extravillous trophoblasts                     Fibroblasts 
##                           20090                           20090 
##   Gastric mucus-secreting cells     Glandular and luminal cells 
##                           20090                           20090 
##                    granulocytes                 Granulosa cells 
##                           20090                           20090 
##                     Hepatocytes                  Hofbauer cells 
##                           20090                           20090 
##                Horizontal cells              Inhibitory neurons 
##                           20090                           20090 
##         Intestinal goblet cells                       Ionocytes 
##                           20090                           20090 
##                   Kupffer cells                Langerhans cells 
##                           20090                           20090 
##                 Late spermatids                    Leydig cells 
##                           20090                           20090 
##                     Macrophages                     Melanocytes 
##                           20090                           20090 
##                Microglial cells                       monocytes 
##                           20090                           20090 
##           Mucus glandular cells               Muller glia cells 
##                           20090                           20090 
##                        NK-cells Oligodendrocyte precursor cells 
##                           20090                           20090 
##                Oligodendrocytes      Pancreatic endocrine cells 
##                           20090                           20090 
##                    Paneth cells               Peritubular cells 
##                           20090                           20090 
##                    Plasma cells       Prostatic glandular cells 
##                           20090                           20090 
##            Proximal enterocytes          Proximal tubular cells 
##                           20090                           20090 
##      Respiratory ciliated cells         Rod photoreceptor cells 
##                           20090                           20090 
##             Salivary duct cells                   Schwann cells 
##                           20090                           20090 
##          Serous glandular cells                   Sertoli cells 
##                           20090                           20090 
##               Skeletal myocytes             Smooth muscle cells 
##                           20090                           20090 
##                   Spermatocytes                   Spermatogonia 
##                           20090                           20090 
##       Squamous epithelial cells        Suprabasal keratinocytes 
##                           20090                           20090 
##            Syncytiotrophoblasts                         T-cells 
##                           20090                           20090 
##                     Theca cells         Thymic epithelial cells 
##                           20090                           20090 
##          Undifferentiated cells 
##                           20090
ct_immune = fread("../Annotation/rna_immune_cell_monaco.tsv.gz")
dim(ct_immune)
## [1] 602700      5
ct_immune[1:5,]
##               Gene Gene name                Immune cell TPM pTPM
## 1: ENSG00000000003    TSPAN6                   basophil  NA  1.2
## 2: ENSG00000000003    TSPAN6  Central memory CD8 T-cell  NA  1.7
## 3: ENSG00000000003    TSPAN6         classical monocyte  NA  0.2
## 4: ENSG00000000003    TSPAN6 Effector memory CD8 T-cell  NA  0.7
## 5: ENSG00000000003    TSPAN6    Exhausted memory B-cell  NA  0.7
summary(ct_immune$TPM)
##    Mode    NA's 
## logical  602700
summary(ct_immune$pTPM)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.00     0.00     3.10    49.74    27.20 96572.50
summary(ct_immune$pTPM[ct_immune$pTPM > 0])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.10     1.70    11.60    67.96    42.80 96572.50
length(unique(ct_immune$`Immune cell`))
## [1] 30
table(ct_immune$`Immune cell`)
## 
##                            basophil           Central memory CD8 T-cell 
##                               20090                               20090 
##                  classical monocyte          Effector memory CD8 T-cell 
##                               20090                               20090 
##             Exhausted memory B-cell               intermediate monocyte 
##                               20090                               20090 
##                         MAIT T-cell               Memory CD4 T-cell TFH 
##                               20090                               20090 
##               Memory CD4 T-cell Th1          Memory CD4 T-cell Th1/Th17 
##                               20090                               20090 
##              Memory CD4 T-cell Th17               Memory CD4 T-cell Th2 
##                               20090                               20090 
##                          myeloid DC                        naive B-cell 
##                               20090                               20090 
##                    naive CD4 T-cell                    naive CD8 T-cell 
##                               20090                               20090 
##                          neutrophil                             NK-cell 
##                               20090                               20090 
##              non-classical monocyte          Non-switched memory B-cell 
##                               20090                               20090 
##                       Non-Vd2 gdTCR                         Plasmablast 
##                               20090                               20090 
##                     plasmacytoid DC                     Progenitor cell 
##                               20090                               20090 
##              Switched memory B-cell                               T-reg 
##                               20090                               20090 
## Terminal effector memory CD4 T-cell Terminal effector memory CD8 T-cell 
##                               20090                               20090 
##                          total PBMC                           Vd2 gdTCR 
##                               20090                               20090
lineage = fread("../Annotation/rna_immune_cell_monaco_cell_types.tsv")
dim(lineage)
## [1] 30  2
lineage[1:2,]
##     Cell_type      Lineage
## 1:   Basophil Granulocytes
## 2: Neutrophil Granulocytes

Check gene expression for each gene set

dim(gene_info)
## [1] 1649    3
for(k in 1:length(gene_sets)){
  if(length(gene_sets[[k]]) < 10) { next }
  
  print(k)
  set_k = paste0("set_", k)
  print(gene_sets[[k]])
  
  genes = gene_info[sym_in_data %in% gene_sets[[k]], gene_symbol]

  n_genes = sum(genes %in% ct_immune$`Gene name`)
  print(sprintf("found %d genes.", n_genes))

  if(n_genes == 0) { next }

  df = ct_immune[`Gene name` %in% genes,]
  dim(df)
  df[1:2,]
  
  stopifnot(all(str_to_lower(df$`Immune cell`) %in% 
                  str_to_lower(lineage$Cell_type)))
  
  mat1 = match(str_to_lower(df$`Immune cell`), 
               str_to_lower(lineage$Cell_type))
  
  df = cbind(df, lineage[mat1,])
  df[1:2,]
  
  df$Cell_type = factor(df$Cell_type, levels = lineage$Cell_type)
  df = df[df$Lineage != "Total PBMC",]
  df$Lineage   = factor(df$Lineage, 
                        levels = setdiff(lineage$Lineage, "Total PBMC"))
  
  p1 = ggplot(df, aes(x=Cell_type, y=log10(pTPM + 0.1), fill=Lineage)) + 
    geom_boxplot() + xlab("Cell type") + ggtitle(set_k) + 
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
    scale_fill_brewer(palette="RdBu")

  print(p1)
}
## [1] 1
##  [1] "AC087623.3" "AMD1"       "CITED2"     "COG5"       "GLTP"      
##  [6] "IFRD1"      "ITGAE"      "MAPRE2"     "RSRP1"      "THAP9-AS1" 
## [11] "WSB1"       "ZFP36L1"    "NEAT1"      "RASGRP1"   
## [1] "found 11 genes."

## [1] 2
##  [1] "CCL4L2"  "ARHGEF3" "ARL4C"   "CLEC16A" "COL6A2"  "COL6A3"  "DDIT4"  
##  [8] "ISG20"   "ITM2A"   "LAIR2"   "LPIN1"   "REXO1"   "TARSL2"  "TRAV4"  
## [15] "XCL1"   
## [1] "found 15 genes."

## [1] 3
##  [1] "ATP5F1A"      "C10orf91"     "COQ8A"        "EPB41L4A-AS1" "SNHG12"      
##  [6] "BROX"         "DMTF1"        "MTERF2"       "PCED1B"       "PPP4R3B"     
## [11] "PTGDS"        "RNF19A"       "SLF2"         "TMEM127"      "TRBV7-6"     
## [16] "TSPAN32"      "VPS13C"      
## [1] "found 14 genes."

## [1] 4
##  [1] "AC245297.3" "AF213884.3" "CXorf40A"   "KCNQ1OT1"   "KMT2E-AS1" 
##  [6] "MATR3-1"    "NPIPB4"     "TRAV8-4"    "TRBV28"     "TRBV7-9"   
## [11] "TRGV7"      "GK5"        "MIAT"       "PARP15"     "SPATA13"   
## [1] "found 8 genes."

## [1] 5
##  [1] "AC044849.1"  "MZF1-AS1"    "TRBV6-1"     "ABCC10"      "AP005482.1" 
##  [6] "CPPED1"      "DENND4B"     "ELMOD3"      "LINC02256"   "MINDY2"     
## [11] "NBEAL2"      "SEC14L1"     "THUMPD3-AS1" "TRBV2"       "TRGV10"     
## [16] "Z93930.2"    "ZNF808"     
## [1] "found 11 genes."

## [1] 6
##  [1] "C12orf45" "CD28"     "GSTM1"    "NR1D1"    "PDE7A"    "TSPYL4"  
##  [7] "ANKRD12"  "ANKRD36C" "NRDC"     "PCSK7"    "PPM1K"    "SLFN12L" 
## [13] "TRIM38"   "UNC13D"   "VTI1A"   
## [1] "found 14 genes."

## [1] 7
##  [1] "CAMK4"  "KCTD7"  "KLRK1"  "MAP3K2" "MZT2B"  "NUAK2"  "WARS2"  "CHD9"  
##  [9] "ERICH1" "IRF9"   "MYO1F"  "NLRC3"  "RASA3"  "STK10" 
## [1] "found 14 genes."

## [1] 8
##  [1] "ALKBH7"  "C6orf48" "EPS8"    "RCAN3"   "TXK"     "GNPTAB"  "KLF2"   
##  [8] "MAN2C1"  "PDE4B"   "RIC3"    "RNF125"  "SCRN3"   "TOB1"    "TRGV4"  
## [15] "ZBP1"   
## [1] "found 14 genes."

## [1] 9
##  [1] "SNHG9"      "ADCY7"      "CCDC84"     "ERBIN"      "FAM78A"    
##  [6] "HPS4"       "KIAA1109"   "NUTM2B-AS1" "RAB27B"     "SIDT1"     
## [11] "SLCO3A1"    "SYTL3"      "TRDV1"      "TTC38"      "ZNF652"    
## [1] "found 12 genes."

## [1] 10
##  [1] "ABCA5"      "AC092683.1" "DDX3Y"      "EIF1AY"     "HECTD4"    
##  [6] "KDM5D"      "RPS4Y1"     "SBNO2"      "TRAV17"     "TTTY15"    
## [11] "UTY"       
## [1] "found 9 genes."

## [1] 11
##  [1] "BEX4"     "CCR7"     "CD40LG"   "CREBL2"   "CRTAM"    "FAM102A" 
##  [7] "FXYD7"    "HIST1H3H" "SDR42E2"  "WDR86"    "EIF4E3"   "PIK3CD"  
## [13] "PYROXD1"  "RUBCN"    "SAMD9L"   "TTC16"   
## [1] "found 15 genes."

## [1] 12
##  [1] "LEF1"    "SELL"    "ZFAND1"  "ACAP3"   "ATAD2"   "CARD11"  "CHD1"   
##  [8] "FAM53B"  "MSI2"    "PIP4K2B" "PKD1"    "RNF157"  "VCAN"   
## [1] "found 13 genes."

## [1] 13
##  [1] "CHRM3-AS2" "RETREG1"   "BMT2"      "FAM133B"   "HRH2"      "IGKV3-15" 
##  [7] "LINC02384" "PUM3"      "S100A12"   "TENT5C"    "THAP5"     "TRAV12-3" 
## [13] "TRAV19"    "TRBV11-2"  "TRBV4-2"  
## [1] "found 12 genes."

## [1] 14
##  [1] "AL118516.1" "APMAP"      "CMC1"       "EFCAB2"     "FAM173A"   
##  [6] "KIR3DL2"    "MATK"       "PITPNC1"    "PRR7"       "IQGAP2"    
## [11] "RLF"        "TBC1D2B"    "TRGC2"      "TUT7"       "XIST"      
## [1] "found 13 genes."

## [1] 15
##  [1] "LINC00402" "SESN2"     "YPEL3"     "CX3CR1"    "GCA"       "LSS"      
##  [7] "MKL1"      "NLRC5"     "ODF3B"     "PATL2"     "PDZD4"     "PEX26"    
## [13] "SZT2"      "ZMIZ2"     "ZNF557"   
## [1] "found 14 genes."

## [1] 16
##  [1] "GALNT11"  "MYLIP"    "ANKRD36"  "ARHGAP10" "ASCL2"    "CD46"    
##  [7] "GON4L"    "MBD5"     "PDE4D"    "PHF14"    "PLAC8"    "RAP1GAP2"
## [13] "SENP7"    "TTC17"   
## [1] "found 14 genes."

## [1] 17
##  [1] "ASL"        "FAM213B"    "AC020659.1" "C2CD3"      "CCDC112"   
##  [6] "DDX60L"     "EPSTI1"     "ETFDH"      "IFI44L"     "IRAK4"     
## [11] "MX2"        "RREB1"      "TRDC"       "XAF1"      
## [1] "found 13 genes."

## [1] 18
##  [1] "LRRN3"     "PAPSS1"    "PDCD4-AS1" "TBCCD1"    "TESPA1"    "TMEM204"  
##  [7] "TRAV21"    "DOCK10"    "ERAP2"     "FGFBP2"    "GALNT3"    "IFI27"    
## [13] "MIGA1"     "OAS2"      "ZNF683"   
## [1] "found 14 genes."

## [1] 19
##  [1] "CLDND1"  "CPNE1"   "GZMK"    "HIKESHI" "KIF9"    "PIK3IP1" "SLC38A2"
##  [8] "UIMC1"   "APOL6"   "TUT4"    "VPS13B"  "ZBTB40"  "ZDHHC5" 
## [1] "found 13 genes."

## [1] 20
##  [1] "NXT2"    "TCEA3"   "ZNF575"  "ABR"     "ARAP1"   "BCL9L"   "C5orf24"
##  [8] "FAM13B"  "GPR141"  "GRK2"    "INPP4A"  "INPP5D"  "SSBP3"   "TRAV27" 
## [15] "ZNF292"  "ZNF493" 
## [1] "found 16 genes."

## [1] 21
##  [1] "AC245014.3" "AL138963.3" "DTHD1"      "LST1"       "NBPF14"    
##  [6] "RGL4"       "TC2N"       "TRAV13-1"   "TRAV8-2"    "TRBV9"     
## [11] "C16orf72"   "GABPB2"     "HECA"       "PCNX1"      "ZNF236"    
## [1] "found 12 genes."

## [1] 22
##  [1] "AC008555.5" "ARRDC2"     "LRRC23"     "NOCT"       "NUP58"     
##  [6] "PGGHG"      "TRAV14DV4"  "TRGV8"      "ZNF749"     "AC116407.2"
## [11] "BICRAL"     "GPR132"     "MYBL1"      "POLR2J3-1"  "TRANK1"    
## [1] "found 12 genes."

## [1] 23
##  [1] "AIF1"    "CD27"    "GIMAP1"  "MSC"     "PCMTD2"  "PDE3B"   "PLK2"   
##  [8] "SLC38A1" "TMEM107" "CASP10"  "COX19"   "DUS1L"   "ITGAM"   "YPEL1"  
## [1] "found 14 genes."

## [1] 24
##  [1] "AC119396.1" "ANXA2R"     "BNIP3L"     "C12orf10"   "CMTM7"     
##  [6] "DYRK4"      "EPHX2"      "HIBADH"     "NT5C3B"     "PPP1R15B"  
## [11] "RGCC"       "SNHG15"     "SNHG8"      "ZNF10"     
## [1] "found 11 genes."

## [1] 25
##  [1] "BTG1"     "HIST1H4C" "LTB"      "RGS10"    "TRBC1"    "CD38"    
##  [7] "GBP1"     "GBP5"     "IFITM2"   "LAG3"     "MIDN"     "NKG7"    
## [13] "SLA2"     "STK17B"  
## [1] "found 14 genes."

## [1] 26
##  [1] "AC004687.1" "CD84"       "CXXC5"      "IER5"       "NCR1"      
##  [6] "RAB33B"     "RPS26"      "TCP11L2"    "TOX"        "TRAV5"     
## [11] "ZFAS1"      "CYTOR"      "MCTP2"      "RUFY2"     
## [1] "found 11 genes."

## [1] 27
##  [1] "IL6R"    "NOSIP"   "ARHGEF9" "CISH"    "CLUH"    "GPHN"    "IL6ST"  
##  [8] "MYO9A"   "NAA25"   "OSM"     "PARP9"   "RALGAPB" "SYNRG"   "ZNF318" 
## [15] "ZNF407" 
## [1] "found 15 genes."

## [1] 28
##  [1] "BTG2"     "CCNB1IP1" "GCSAM"    "LDLRAP1"  "SERINC5"  "TCF7"    
##  [7] "WDR54"    "ABCA7"    "DIAPH2"   "EHMT1"    "HIPK3"    "KIAA2026"
## [13] "NEK9"     "USP16"   
## [1] "found 14 genes."

## [1] 29
##  [1] "COA1"      "SLC2A3"    "TMIGD2"    "AP3B1"     "AP3M2"     "CEMIP2"   
##  [7] "EHBP1L1"   "FGL2"      "HEXDC"     "INO80D"    "LINC02446" "MAF"      
## [13] "OXNAD1"    "TMEM131L" 
## [1] "found 13 genes."

## [1] 30
##  [1] "AC004854.2" "AC087239.1" "AL136454.1" "AL627171.1" "BX284668.6"
##  [6] "LINC02273"  "MCUB"       "MFNG"       "RGS1"       "Z93241.1"  
## [11] "CRYBG1"     "GCN1"       "GDPD5"      "MARF1"      "MICAL2"    
## [16] "PSMA3-AS1"  "ZNF83"     
## [1] "found 9 genes."

## [1] 31
##  [1] "AC012645.3" "AC016405.3" "AC020911.2" "AC027644.3" "AC083798.2"
##  [6] "AC091271.1" "AK5"        "AL135791.1" "ARF4-AS1"   "LINC01465" 
## [11] "PRAG1"      "TRAV12-2"   "TRAV8-3"    "TRBV3-1"    "TRBV6-2"   
## [16] "ZNF862"    
## [1] "found 7 genes."

## [1] 32
##  [1] "CARHSP1"  "CRLF3"    "CYB561A3" "NT5DC1"   "ARAP2"    "ARHGAP30"
##  [7] "DGKD"     "HIPK1"    "HSH2D"    "NCKAP1L"  "NECAP1"   "PRR14L"  
## [13] "R3HCC1L"  "SETX"     "TMEM8A"  
## [1] "found 15 genes."

## [1] 33
##  [1] "AC015982.1" "AC025164.1" "AC083880.1" "AL121944.1" "AL451085.1"
##  [6] "ARMH1"      "ATP2B1-AS1" "CSKMT"      "CXorf57"    "HIPK1-AS1" 
## [11] "ILF3-DT"    "INPP4B"     "OSER1-DT"   "TRABD2A"   
## [1] "found 5 genes."

## [1] 34
##  [1] "SESN1"       "TRBV20-1"    "TRBV7-2"     "ADGRE5"      "AKNA"       
##  [6] "ARHGAP45"    "ENOSF1"      "IQCG"        "NFATC3"      "OGA"        
## [11] "PIK3R5"      "SLC16A1-AS1" "SUSD6"       "TRAV9-2"     "XCL2"       
## [1] "found 14 genes."

## [1] 35
##  [1] "AOAH"        "FAM118A"     "HLA-DMB"     "KLRC3"       "KLRF1"      
##  [6] "MTRNR2L8"    "TRAV38-2DV8" "TRBV6-5"     "ADGRG1"      "FAM126B"    
## [11] "FCRL6"       "LILRB1"      "TMEM181"    
## [1] "found 12 genes."

## [1] 36
##  [1] "AC007952.4" "AC103591.3" "AL357060.1" "BBS9"       "C6orf99"   
##  [6] "HELQ"       "IGKV3-20"   "IGLV1-44"   "INTS6L"     "JAML"      
## [11] "LINC00649"  "TNFRSF25"   "TRAV1-2"    "TRAV3"      "TRGV5"     
## [1] "found 10 genes."

## [1] 37
##  [1] "AC007384.1" "AC025171.3" "AL139246.5" "ARL4A"      "GADD45B"   
##  [6] "ID1"        "NR4A2"      "NR4A3"      "AC016831.7" "ELMO1"     
## [11] "FAM160B1"   "GPRIN3"     "NARFL"      "PRR5L"      "ZBTB20"    
## [1] "found 10 genes."

## [1] 39
##  [1] "ARHGAP9" "EOMES"   "FCRL3"   "TRG-AS1" "TRGV9"   "ATAD2B"  "CCL4"   
##  [8] "CREBZF"  "CST7"    "CTSW"    "GPR174"  "ITGAL"   "MYO1G"   "NFE2L1" 
## [1] "found 13 genes."

## [1] 40
##  [1] "ADAM17"  "ADHFE1"  "CD226"   "CHD6"    "ETNK1"   "LONP2"   "MPPE1"  
##  [8] "SEMA4D"  "ST8SIA4" "SUSD1"   "TGFBR3"  "WDTC1"  
## [1] "found 12 genes."

Session information

gc()
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  4174450 223.0   11673970 623.5         NA 11673970 623.5
## Vcells 20587554 157.1   61723001 471.0      65536 61693859 470.7
sessionInfo()
## R version 4.2.3 (2023-03-15)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] stringr_1.5.0     limma_3.54.2      tidyr_1.3.0       ggpubr_0.6.0     
## [5] ggplot2_3.4.2     data.table_1.14.8
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.10            png_0.1-8              Biostrings_2.66.0     
##  [4] digest_0.6.31          utf8_1.2.3             R6_2.5.1              
##  [7] GenomeInfoDb_1.34.9    backports_1.4.1        stats4_4.2.3          
## [10] RSQLite_2.3.1          evaluate_0.20          httr_1.4.6            
## [13] pillar_1.9.0           zlibbioc_1.44.0        rlang_1.1.0           
## [16] rstudioapi_0.14        car_3.1-2              jquerylib_0.1.4       
## [19] blob_1.2.4             R.oo_1.25.0            R.utils_2.12.2        
## [22] S4Vectors_0.36.2       rmarkdown_2.21         labeling_0.4.2        
## [25] RCurl_1.98-1.12        bit_4.0.5              munsell_0.5.0         
## [28] broom_1.0.4            compiler_4.2.3         xfun_0.39             
## [31] pkgconfig_2.0.3        BiocGenerics_0.44.0    htmltools_0.5.5       
## [34] tidyselect_1.2.0       KEGGREST_1.38.0        GenomeInfoDbData_1.2.9
## [37] tibble_3.2.1           IRanges_2.32.0         fansi_1.0.4           
## [40] crayon_1.5.2           dplyr_1.1.2            withr_2.5.0           
## [43] R.methodsS3_1.8.2      bitops_1.0-7           grid_4.2.3            
## [46] jsonlite_1.8.4         gtable_0.3.3           lifecycle_1.0.3       
## [49] DBI_1.1.3              magrittr_2.0.3         scales_1.2.1          
## [52] cli_3.6.1              stringi_1.7.12         cachem_1.0.7          
## [55] carData_3.0-5          farver_2.1.1           XVector_0.38.0        
## [58] ggsignif_0.6.4         bslib_0.4.2            generics_0.1.3        
## [61] vctrs_0.6.2            RColorBrewer_1.1-3     org.Hs.eg.db_3.16.0   
## [64] tools_4.2.3            bit64_4.0.5            Biobase_2.58.0        
## [67] glue_1.6.2             purrr_1.0.1            abind_1.4-5           
## [70] fastmap_1.1.1          yaml_2.3.7             AnnotationDbi_1.60.2  
## [73] colorspace_2.1-0       rstatix_0.7.2          memoise_2.0.1         
## [76] knitr_1.44             sass_0.4.5